#### PRODUCES SIMULATIONS FOR WHAT MIGHT HAVE GENERATED THE OBSERVED BVOTING PATTERNS
#### GIVEN THE ESTIMATES OF INDIVIDUAL EFFECTS WE HAVE
#### PRODCUES GRAPHS USED IN BH PRESENTATION AND BPSR PIECE

#### USes actual extracts of actual 2010 and 2014 data on coverage as starting point


####### Thinking about what happened *before*###########
###### Possible comparisons between 2010 and 2014 ######

#### Make predictions for 2014
load("data_replicsimulation2014.RData")
s2014$coverage2014 <- s2014$alcance.2014
s2014$poverty2014 <- pmax(s2014$targetalcance.2014,s2014$coverage)
### Compute commulative vote shares ###
s2014 <- s2014[order(s2014$coverage2014,decreasing=F),]
s2014$cumvote <- cumsum(s2014$dilma.2014)
s2014$cumvs <- s2014$cumvote/sum(s2014$dilma.2014,na.rm=T) 
### Compute reverse commulative vote shares ###
s2014 <- s2014[order(s2014$coverage2014,decreasing=T),]
s2014$revcumvote <- cumsum(s2014$dilma.2014)
s2014$revcumvs <- s2014$revcumvote/sum(s2014$dilma.2014,na.rm=T)

### Predict votes and predicted commulative votes/vs for 2014 #####
elbow <- .35
s2014$poor.nb <- .075+ifelse(s2014$coverage2014<elbow
									,0.50164583-.11
									,0.50164583-.11-(elbow*.46)+s2014$coverage2014*.46)
s2014$poor.b <- 0.12342628+ s2014$poor.nb
s2014$non.poor <- #.1+0.50164583-.15-(.35*.5)+s2014$coverage2014*.5
					.23+s2014$coverage*.52
s2014$votes.poor <- s2014$coverage*s2014$poor.b + (s2014$poverty2014-s2014$coverage2014)*s2014$poor.nb
s2014$votes.non.poor <- (1-s2014$poverty)*s2014$non.poor
s2014$predvs <- s2014$votes.poor + s2014$votes.non.poor
s2014$predvotes <- s2014$predvs*s2014$valid.2014
s2014 <- s2014[order(s2014$coverage2014,decreasing=F),] #cummulative
s2014$predcumvote <- cumsum(s2014$predvotes)
s2014$predcumvs <- s2014$predcumvote/sum(s2014$predvotes,na.rm=T) 


### Load data and make predictions for 2010 ###
load("data_replicsimulation2010.RData")
s$coverage2010 <- s$alcance.2010
s$poverty2010 <- pmax(s$targetalcance.2006,s$coverage2010) 
### Compute commulative vote shares ###
s <- s[order(s$coverage2010,decreasing=F),]
s$cumvote <- cumsum(s$dilma.2010)
s$cumvs <- s$cumvote/sum(s$dilma.2010,na.rm=T)
### Compute reverse commulative vote shares ###
s <- s[order(s$coverage2010,decreasing=T),]
s$revcumvote <- cumsum(s$dilma.2010)
s$revcumvs <- s$revcumvote/sum(s$dilma.2010,na.rm=T)
### Predict votes and predicted commulative votes/vs for 2010 #####
elbow <- .35
s$poor.nb <- .1+ifelse(s$coverage2010<elbow 
									,0.54889952-.11
									,0.54889952-.11-(elbow *.46)+s$coverage2010*.46)
s$poor.b <- 0.09755105+ s$poor.nb
s$non.poor <- #.1+0.54889952-.15-(.35*.5)+s$coverage2010*.5
			.33+s$coverage*.35#+s$coverage^2*.1
s$votes.poor <- s$coverage2010*s$poor.b + (s$poverty2010-s$coverage2010)*s$poor.nb
s$votes.non.poor <- (1-s$poverty2010)*s$non.poor
s$predvs <- s$votes.poor + s$votes.non.poor
s$predvotes <- s$predvs*s$valid.2010
s <- s[order(s$coverage2010,decreasing=F),] #cummulative
s$predcumvote <- cumsum(s$predvotes)
s$predcumvs <- s$predcumvote/sum(s$predvotes,na.rm=T) 

#### Merge 2010 and 2014
tmp <- merge(s2014,s,by="codeibge",all=T,suffixes=c("2014","2010"))

#### Plots for both years ######
#### Vote share #####
par(mfrow=c(2,2))
tmp <- tmp[order(tmp$coverage2010),]
plot(tmp$coverage2010,tmp$dilma.vs.2010,type="n",ylim=c(0,1),main="2010")
lines(tmp$coverage2010,tmp$poor.nb2010,col=gray(0.5))
lines(tmp$coverage2010,tmp$poor.b2010,col=gray(0.5))
lines(tmp$coverage2010,tmp$non.poor2010,col=gray(0.5))
#lines(lowess(na.omit(cbind(tmp$coverage2010,tmp$non.poor2010))),lty=2,col=gray(0.5))
lines(lowess(na.omit(cbind(tmp$coverage2010,tmp$dilma.vs.2010))))
lines(lowess(na.omit(cbind(tmp$coverage2010,tmp$predvs2010))),col=2,lty=2)

tmp <- tmp[order(tmp$coverage2014),]
plot(tmp$coverage2014,tmp$dilma.vs.2014,type="n",ylim=c(0,1),main="2014")
lines(tmp$coverage2014,tmp$poor.nb2014,col=gray(0.5))
lines(tmp$coverage2014,tmp$poor.b2014,col=gray(0.5))
lines(lowess(na.omit(cbind(tmp$coverage2014,tmp$non.poor2014))),lty=2,col=gray(0.5))
lines(lowess(na.omit(cbind(tmp$coverage2014,tmp$dilma.vs.2014))))
lines(lowess(na.omit(cbind(tmp$coverage2014,tmp$predvs2014))),col=2,lty=2)

#### Ccummulative vote share ####
plot(tmp$coverage2010,tmp$cumvote2010/1000000,pch=".",xlim=c(0,1)
		,ylab="Millions of Votes",xlab="BF Coverage",main="2014",col=1)
points(tmp$coverage2010,tmp$predcumvote2010/1000000,pch=".",col=2)
plot(tmp$coverage2010,tmp$cumvote2010/1000000,type="n",xlim=c(0,1)
		,ylab="Millions of Votes",xlab="BF Coverage",main="2014")#to force same scale
points(tmp$coverage2014,tmp$cumvote2014/1000000,pch=".",col=1)
points(tmp$coverage2014,tmp$predcumvote2014/1000000,pch=".",col=2)


#### Variations in plots for BPSR paper ####
jpeg(file="fig-BPSRsimulations20102014.jpg"
		,width=10,height=5.4,units="in",type="quartz",res=300)
par(mfrow=c(1,2),mar=c(4,4,3,1))
tmp <- tmp[order(tmp$coverage2010),]
plot(tmp$coverage2010,tmp$dilma.vs.2010,type="n",ylim=c(0,1),main="2010",ylab="",xlab="")
mtext(side=1,line=2.5,"Bolsa Família Coverage")
mtext(side=2,line=2.5,"Prob of Voting for Dilma")
lines(tmp$coverage2010,tmp$poor.nb2010,col=gray(0.5))
lines(tmp$coverage2010,tmp$poor.b2010,col=gray(0.5))
lines(tmp$coverage2010,tmp$non.poor2010,col=gray(0.5))
N <- which.min(tmp$coverage2010)
text(x=rep(-.05,2),y=s[N,c("poor.nb","poor.b")]+c(-.025,+.017)
		,labels=c("Poor Non-BF","BF"),pos=4,xpd=NA,cex=.9)
text(.05,s$non.poor[N],labels="Non-Poor",pos=4,cex=.9)
arrows(x0=tmp$coverage2010[N],x1=tmp$coverage2010[N]
		,y0=tmp$poor.nb2010[N],y1=tmp$poor.b2010[N],code=3,length=.1)
text(s$coverage[N],s$poor.nb[N]+.049,labels="=0.098",pos=4,cex=.9)

tmp <- tmp[order(tmp$coverage2014),]
plot(tmp$coverage2014,tmp$dilma.vs.2014,type="n",ylim=c(0,1),main="2014",ylab="",xlab="")
mtext(side=1,line=2.5,"Bolsa Família Coverage")
mtext(side=2,line=2.5,"Prob of Voting for Dilma")
lines(tmp$coverage2014,tmp$poor.nb2014,col=gray(0.5))
lines(tmp$coverage2014,tmp$poor.b2014,col=gray(0.5))
lines(tmp$coverage2014,tmp$non.poor2014,col=gray(0.5))
N <- which.min(tmp$coverage2014)
text(x=rep(-.05,2),y=tmp[N,c("poor.nb2014","poor.b2014")]+c(-.025,+.017)
		,labels=c("Poor Non-BF","BF"),pos=4,xpd=NA,cex=.9)
text(.05,tmp$non.poor2014[N],labels="Non-Poor",pos=4,cex=.9)
arrows(x0=tmp$coverage2014[N],x1=tmp$coverage2014[N]
		,y0=tmp$poor.nb2014[N],y1=tmp$poor.b2014[N],code=3,length=.1)
text(tmp$coverage2014[N],tmp$poor.nb2014[N]+.06,labels="=0.123",pos=4,cex=.9)
dev.off()




plot(tmp$coverage2010,tmp$cumvote2010/1000000,pch=".",xlim=c(0,1),ylab="",xlab="")
mtext(side=1,line=2.5,"BF Coverage")
mtext(side=2,line=2.5,"Millions of Votes")
points(tmp$coverage2014,tmp$cumvote2014/1000000,pch=".",col=1)
text(.6,30,labels="2014")
text(.56,40,pos=2,labels="2010")

plot(tmp$coverage2010,tmp$dilma.vs.2010,type="n",ylim=c(0,1),ylab="",xlab="")
lines(lowess(na.omit(cbind(tmp$coverage2010,tmp$dilma.vs.2010))))
lines(lowess(na.omit(cbind(tmp$coverage2014,tmp$dilma.vs.2014))))
mtext(side=1,line=2.5,"BF Coverage")
mtext(side=2,line=2.5,"Average Vote Share")
text(.6,.6,labels="2014")
text(.56,.69,pos=2,labels="2010")



